Optimally coherent sets in geophysical flows: A new approach to delimiting the 
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The "edge" of the Antarctic polar vortex is known to behave as a barrier to the meridional (pole- 
ward) transport of ozone during the austral winter. This chemical isolation of the polar vortex 
from the middle and low latitudes produces an ozone minimum in the vortex region, intensifying 
the ozone hole relative to that which would be produced by photochemical processes alone. Ob- 
servational determination of the vortex edge remains an active field of research. In this letter, 
we obtain objective estimates of the structure of the polar vortex by introducing a new technique 
based on transfer operators that aims to find regions with minimal external transport. Applying 
this new technique to European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-40 
three-dimensional velocity data we produce an improved three-dimensional estimate of the vortex 
location in the upper stratosphere where the vortex is most pronounced. This novel computational 
approach has wide potential application in detecting and analysing mixing structures in a variety 
of atmospheric, oceanographic, and general fiuid dynamical settings. 



I. INTRODUCTION 

Stratospheric ozone in the Southern Hemisphere high 
latitudes has decreased dramatically since the early 
1970s. This long-term trend has been attributed to a 
combination of natural and anthropogenic factors [H-IJ]- 
In particular, it has been discovered that the ozone de- 
pletion in the lower stratosphere of the Southern Hemi- 
sphere is particularly pronounced, due in part to a strong 
barrier to meridional transport between middle and high 
latitudes during the austral winter and early spring 
Barriers such as these, which often coexist with turbulent 
mixing, play a major role in the dynamics of the strato- 
sphere. The polar vortex is a known strong barrier to 
transport, enclosing a persistent, non-dispersive, coher- 
ent region over the high latitudes. Our aim in this letter 
is to precisely determine the spatial location and move- 
ments of this coherent region, improving significantly 
over existing methods of estimation. Our study focuses 
on the upper stratosphere where the polar vortex is best 
developed. 

It is common meteorological practice to diagnose the 
polar vortex edge at the position of maximum meridional 
gradient of potential vorticity (PV). Potential vorticity is 
a quantity combining measures of circulation and strat- 
ification which is materially conserved for adiabatic, in- 
viscid flow (both of which are good approximations in 
stratospheric flow over timescales of a week or two). It 
can be shown that strong PV gradients produce a "restor- 
ing force" inhibiting meridional motion of air parcels Q . 
Nevertheless, PV gradients alone provide only indirect 
measures of mixing barriers. In contrast, the present 
study characterises regions of minimal mixing directly in 
terms of the transport properties of the observed strato- 
spheric flow. We present an innovative new mathematical 
technique to determine the polar vortex location at dif- 
ferent times, directly as coherent structures in observed 



velocity fields. Lagrangian PV-based measures of the 
vortex such as those presented in [l^ are complicated 
by the fact that PV is generally a noisy field (as vorticity 
is the curl of the velocity field). The velocity field is gen- 
erally much smoother; barriers to mixing estimated from 
the velocity field can be expected to be less sensitive to 
(poorly-observed) small-scale features of the flow. 

Our new mathematical approach for detecting mini- 
mal transport structures with high accuracy has a broad 
range of potential applications to geophysical fluids. For 
example, transport properties in other long-lived atmo- 
spheric coherent structures such as blocking highs are of 
interest. There is also increasing interest in the transport 
properties of mesoscale (on the order of 10 to 100 km in 
diameter) ocean eddies and their influence on biological 
processes within the upper, sunlit part of the water col- 
umn 0,0]. 



II. INPUT DATA AND NON-AUTONOMOUS 
FLOW 

Our input data consists of three-dimensional veloc- 
ity flelds obtained from the ECMWF ERA-40 data set 
( |http://data.ecmwf.int/data/index.html| . The data is 
on a three-dimensional grid with 2.5 degree resolution 
in the latitude and longitude direction (144 by 73 grid 
points over the Southern Hemisphere). Vertical coordi- 
nates are in units of hPa, with data provided at 5 pressure 
levels (3, 5, 7, 10, 20 hPa). We use 62 days of 6-hourly 
velocity flelds from August 1 to September 31 in 1999. 
The velocity fields will be interpolated linearly in space 
and in time; thus we can only aim to detect features at 
the resolution of the data provided. While we recognize 
that there may be biases in the reanalysis data, particu- 
larly in the upper pressure levels near the model's upper 
boundary, the purpose of this study is to demonstrate the 
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ability of the transfer operator approach to characterize 
coherent sets in highly unsteady flows. A climatology of 
the polar vortex would require a more careful considera- 
tion of the dataset under consideration 

Our interest is focused on the Lagrangian dynamics in 
the higher latitudes of Southern Hemisphere. Therefore, 
we will work on the phase space X = x [—90°, —30°] x 
D, where is a circle parameterized from 0° to 360° and 
D = [3, 20] denotes the range of pressure in hPa. The 
Lagrangian motion of passive tracers is represented by 
their trajectories x{t) := $(x,t;T), where the flow map 
$:XxRxR^>Xisa function of time r and gives 
the terminal point $(x,i;r) in X of a particle initially 
located at x G X at time t and flowing for r time units. 
The flow map ^{xq, to; r) is obtained as a solution of the 
nonautonomous ODE ^ = f{x{t),t) with initial condi- 
tion x(io) = ^{xq, to] 0), where f{x{t), t) in this report is 
the prescribed velocity data. 



III. ALMOST-INVARIANT SETS, COHERENT 
PAIRS, AND TRANSFER OPERATORS 

We shall be interested in finding a pair of sets At , At+r 
at times t and t -(- t so that ^{At+r^t + —t) w At- 
Moreover, this pair of sets should retain this property 
even when some diffusion is added to the system. Let /i 
be a probability measure that is preserved by the fiow at 
all times. We call At, At+r a {po,t,T) — coherent pair if 

p^{At,At+r) ■■= ^^iAt n <i>{At+r,t + t; -T))/^i{At) > Pa, 

, (1) 

and p{At) = /i(At+r)- The condition on addition of 
diffusion is crucial. Clearly, there are many {pQ,t,T)- 
coherent pairs according to the above definition without 
diffusion. One may simply select any set At d X and 
define At+r = ^{At,t;T) to produce a (1, T)-coherent 
pair. In chaotic systems, such an image set At+r is likely 
to be significantly less regular than At because of stretch- 
ing and folding. We are seeking (po, i, T)-coherent pairs 
with both sets regular. The requirement that ([T]) hold 
even under diffusion acts as a selection principle, remov- 
ing irregular sets, and selecting pairs that are robust to 
perturbation. At a certain level of diffusion, we may ask 
to find a coherent pair that maximises po, and may ex- 
pect a unique such pair. 

To identify sets satisfying ([J), we use a transfer oper- 
ator Vt,T ■ L^{X,m) O defined by 



g{^{x,t 



-r)) 



I det D^{^{x, t + r; -t), t; r) 



(2) 



where m is the normalized Lebesgue measure on X. In 
particular, if g{x) is a density of passive tracers at time 
t, Vt.rgix) provides their density at time t + t induced 
by the flow 

In the autonomous setting, eigenfunctions / of Vt^r 
(= Vr for all t) corresponding to positive ei gen values 
Awl were used to find almost-invariant sets IsMllll. The 



key point of difference between these prior studies and the 
present work is that the sets studied previously do not 
move significantly over the time duration studied, while 
our present work seeks highly mobile coherent sets, which 
are far from being almost-invariant. The new theory and 
numerics we introduce in the next section are specifically 
designed for nonautonomous or time-dependent systems. 



IV. NUMERICAL APPROACH 

We partition X into a grid of boxes {Bi, . . . , B„}. The 
pressure extents of the boxes are either 3-5, 5-7, 7-10, or 
10-20 hPa. Each pressure layer consists of 6605 boxes 
of approximately equal cross-sectional area in the lati- 
tude/longitude directions, leading to n = 26420 boxes 
in total. To numerically approximate the transfer opera- 
tor Vt^T, we construct a finite-dimensional approximation 
based on Ulam's approach [12]: 



m{B,r\^{Bj,t + T;-T)) 



(3) 



where m is a normalised volume measure in 
(lat, Ion, pressure) coordinates. The entry P('^)(i)ij- 
may be interpreted as the probability that a point 
selected uniformly at random in Bi at time t will 
be in Bj at time t + t. The discretisation naturally 
produces a diffusion at the level of box diameters. As 
our boxes are of approximately the same dimensions as 
the distances between neighbouring ERA-40 data points, 
it is unnecessary to impose additional diffusion. If our 
boxes were significantly smaller the distances between 
neighbouring data points, it is possible that spurious 
fine features below the resolution supported by the data 
could appear; in such a situation, additional diffusion 
would be required to remove spurious fine features. 
We estimate P^^'(Oij by 

pM(i),,, « #{£ : y,,, e B„$(j/,,,,i;T) G B,}/Q, (4) 

where yi^i, i = 1,...,Q are uniformly distributed test 
points in Bi and ^{yi,£,t;T) is obtained via a numerical 
integration. We set Q = 147 in our experiments and 
calculate <^{yi^i,t;T) using the standard Runge-Kutta 
method with step size of 3/4 hours. Linear interpolation 
is used to evaluate the velocity vector of a tracer lying 
between the data grid points in the longitude-latitude- 
pressure coordinate. In the temporal direction the data 
is affinely interpolated independently in the longitude, 
latitude and pressure level directions. The step size of 
3/4 hours is small enough to guarantee that a tracer will 
usually not fiow to a neighboring data grid set; this limits 
the numerical integration error. 

We assume that the mass density of particles in X is 
at equilibrium and denote the fractional mass of particles 
contained in Bi by pi. Thus X^ILi Pi — ^- We construct a 
reverse time transition matrix from time i-f r to t denoted 
P(-)(t) as P(^HO.j = P^"H%.PjM- 
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Introduce a weighted inner product (x, y)p := 
Eti x^y^P^■ One has (xPM(t),y)p = {x,yP^^Ht))p for 
all x,y G M". 

Our new approach to finding a coherent pair is intu- 
itively based upon seeking a solution to 

max (-P'-n;),^P^-HO)._ 

w£{±l}" {w,W)p 

We think of At :— Ui:^.=ii?i and := Ui:wi=-iBi as 
a coherent partition of X. The numerator in ([5]) repre- 
sents the size of the forward image of the vector w. If 
there is little transport from At to $(j4j,i;r) and from 
A^ to $(At,t;r) (so Af,$(Af,t;T) and A^t,^{A^,t;T) 
are both coherent pairs), this numerator will be large. 
To produce non-trivial partitions, we may need to place 
lower bounds on the masses of both At and Af. Such 
a balanced bisection problem is combinatorially hard 
to solve. Therefore we remove the discrete condition 
w e {±1}", allowing w to float freely in M". To ef- 
fect a balancing of mass between positive and negative 
components of w, we insert the condition (w, v)p = 0, for 
some nonnegative test vector v e R" . We will see shortly 
that the correct choice of v is the minimizer of the central 
inner product. Thus, we have 

{wV(^){t),wV(-\t))p 
mm max ; ^ . (o) 

Letting Dij = Sijpt and noting {w, v)p = 
{wD^/'^ ,vD^/'^)2, this is easily solved by computing 
the second largest singular value s oi D~^/^'P^'^\t)D^/^ . 
Denote the corresponding left singular vector by y 
(under multiplication on the right). The maximizing 
w = w{t) is constructed as w (t) = yD-i/z. The 
minimizing v turns out to be uD^^/^ where u is the 
leading left singular vector of D-^^^P'^'^'^ {t)D^/^ . We 
also construct z as the corresponding right singular 
vector and set w'{t + t) = zD^^^^. We assume that 
w{t),'w'{t + t) are normalised so that {w{t),w{t))p = 1 
and {w'{t + t), w'{t + t))p = 1. 
One now has: 

1. u;(i)pM(<) ^ sw'{t + T), 

2. ii;'(i-f T)PM(i) = sw{t), 

3. {w{t)V'^^\t),w{t)Vi^Ht))p = s\ 

Choosing v via the minimization in (|6]) ensures that w 
has the transformation properties 1. and 2. above, which 
are crucial to the definition of coherent sets. 

We now extract a coherent pair At and At+r from 
a pair of vectors w{t) and w'lt + r). We create sets 
that are unions of boxes with w-values above certain 
thresholds. Define A+(c) := Ui:u.(t)>c^' ^^^d A+^^(c) 
y^i■.w'(t+Ty>cB^, c G M. Denote ^„(l+(c)) = Y.i:w(t)>cPi 

and ^ln{At^r{c)) = T.i,w-[t+r)>cPi- 



For At = Uie/t ^"^^ = ^ieit+^ define 

The quantity p„ measures the discretised coherence for 
the pair At^At+r- Our procedure is summarised below: 

1. Let 7y(c) = argminc'eB|Mn(^t'"(c)) - ^ln{Af^^{c'))\. 
This is to enforce /i„(A() = ^niAt+r)- 

2. Set c* = argmaxceBPn(A^(c), At+^^(?7(c))). The 
value of c* is selected to maximize the coherence. 

3. Define At := A+(c*) and At+r :== A+^^{ri{c*)). 

We remark that one has to ensure that the sign of w{t) 
and w'(t-\-T) manifest the same "parity", i.e., the salient 
features of w{t) and w'{t + r) to be extracted must have 
the same sign. It may thus be necessary to multiply one 
of w[t) or w'{t + t) by —1. 

The major computational cost is the construction of 
Pf"^'. The calculation of large singular values and cor- 
responding singular vectors is relatively quick, as pf'^) 
is very sparse and iterative methods for sparse matri- 
ces may be used. The construction of P^"^) requires nu- 
merical integration of Q • n trajectories for a flow dura- 
tion of r time units. The trajectory computations are 
of course highly parallisable, and further computational 
savings might be made by reusing already computed tra- 
jectory segments to link with new trajectories when the 
latter pass nearby. 

V. NUMERICAL RESULTS 

We computed the SVD of D-i/2p(*)(T-)Di/2 at t = 14 
with r = 14 days to obtain the left (resp. right) sin- 
gular vectors y (resp. z) and hence w(14) and w'{2d,). 
Figure [T] illustrates the vectors w(14) and w'(28) with 
the components monotonically rescaled to uniformly dis- 
tributed values between and 1. The highlighted part 
of these vectors describes the most coherent pair of sets. 
We now threshold w(14) and w'(28) using the algorithm 
described above to extract the corresponding pair of co- 
herent sets; see Figure [H We find the optimal coherence 
ratio is /9„(Ai4,A28) ~ 0.7902; this means that about 
21% of the mass in An on August 14, 1999 falls outside 
on August 28, 1999. 

Interestingly, our coherent pair has a "hole" over the 
south pole. Further calculations have revealed that the 
reason for this is that around twice as many particles in 
this vertical hole on August 14 have exited the slice 3-20 
hPa vertically by August 28 when compared to similar 
vertical exits of particles starting in the identified coher- 
ent set on August 14. Thus, this inner part of the vor- 
tex is less coherent and excluded from our coherent pair. 
This hole may be an artifact of the reanalysis data, al- 
though it is consistent with evidence of very strong polar 
descent in this region [T3| . 
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We now compare our coherent pair of sets to sets de- 
fined by contours of potential vorticity (PV) . A common 
approach, developed in [HI is to define the vortex 
boundary as the isoline of the largest gradient of PV 
w.r.t. the equivalent latitude. We employ this approach 
to define potential coherent pairs at t = 14 and t — 28. 
We additionally enforce the constraint that the mass en- 
closed by a PV isocontour at t = 28 is approximately 
equal to the mass of the set enclosed by the determined 
PV isocontour at t = 14. The computational cost of the 
PV approach is NARATIP, PLEASE ADD MATERIAL 
ON COMPUTATIONAL COST. 

The two-dimensional plots of PV-determined coherent 
pairs at t = 14 and t — 28 are compared with the co- 
herent sets in Figure [2] To estimate the transport of 
particles from the inside the set at t = 14 to outside the 
set at i = 28, we use a method similar to the contour 
crossing method introduced in fl^. The tracer particle 
is considered to be outside the boundary if its potential 
vorticity is larger than that of the boundary. Note that 
the contour crossing method is originally developed to es- 
timate the transport on the 2D isentropic surface but we 
would like to extend its utility to estimate the transport 
across the boundary surface. Therefore, we interpolate 
the PV at the final time {t = 28) to obtain the PV at the 
particle's final position. We also interpolate the PV of 
the boundaries of the set at t = 28 along the pressure co- 
ordinate to determine the boundary at the pressure level 
the advected particle resides in. This calculation shows 
that the fraction of particles initially inside the surface at 
t — 14 remains inside the boundary surface at t = 28 is 
approximately 0.7204. Thus our transfer operator based 
approach yields coherent pairs with 9.69% greater co- 
herence. Moreover, our approach is able to detect finer 
structures, including a vertical hole near the south pole. 

VI. CONCLUSIONS 

The Antarctic polar vortex is a well-known feature 
of the austral wintertime stratosphere separating polar 
and midlatitude air masses. The strong barrier to trans- 
port at the vortex edge plays an important role in ozone 
dynamics, particularly the development of the Southern 
Hemisphere ozone hole in austral spring. Diagnosis of the 
vortex edge from observations is a challenging problem 
that remains a subject of active research. 

Previous approaches to this problem have been based 
on kinematic (following the advection of some tracer) or 
dynamic (considering gradients of PV) arguments. We 
presented a new kinematic method of accurately estimat- 
ing the three-dimensional location of the vortex. This 
new method uses the velocity field to diagnose "optimally 
coherent pairs" and was able to determine a significantly 
more accurate estimate of transport barriers, with almost 



10% less external transport from the identified vortex re- 
gion than the PV-based estimate. Future, more detailed 
studies will include an investigation of the climatology 
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FIG. 1: Left column: the vector w{14) shown for the pressure levels 
3-5, 5-7, 7-10 and 10-20 hPa. The 4x6605 = 26420 components 
of w)(14) have been mapped to the values 1/26420,2/26420,. . . ,1, 
preserving their order. Center column: w{28). Right column: Po- 
tential vorticity (K'm?kg~^sec~^) at levels 3, 5, 7 and 10 hPa on 
August 14, 1999 (i = 14). 



of the polar vortex on isentropic surfaces throughout the 
stratosphere. 

Our new computational approach for detecting mini- 
mal transport structures with high accuracy has a broad 
range of potential application to studies of transport and 
mixing in the atmosphere and ocean, and in general fiuid 
dynamics settings. 
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FIG. 2: Comparison between coherent pair and PV surface 
boundary. The optimal coherent set A14 and ^428 obtained from 
thresholding the vectors w{l4) and ui'(28). The coherent ratio 
Pij.{Ai4, A2s) ~ 0.7902. The blue curve in each plot shows the PV 
surface boundary obtained from the maximum PV gradient w.r.t 
equivalent latitude 
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